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ABSTRACT 

The couplings between supermassive black-hole binaries and their environments within galactic nuclei have 
been well studied as part of the search for solutions to the final parsec problem. The scattering of stars by the 
binary or the interaction with a circumbinary disk may efficiently drive the system to sub-parsec separations, 
allowing the binary to enter a regime where the emission of gravitational waves can drive it to merger within 
a Hubble time. However, these interactions can also affect the orbital parameters of the binary. In particular, 
they may drive an increase in binary eccentricity which survives until the system’s gravitational-wave signal 
enters the pulsar-timing array band. Therefore, if we can measure the eccentricity from observed signals, we 
can potentially deduce some of the properties of the binary environment. To this end, we build on previous 
techniques to present a general Bayesian pipeline with which we can detect and estimate the parameters of an 
eccentric supermassive black-hole binary system with pulsar-timing arrays. Additionally, we generalize the 
pulsar-timing array J>-statistic to eccentric systems, and show that both this statistic and the Bayesian pipeline 
are robust when studying circular or arbitrarily eccentric systems. We explore how eccentricity infiuences the 
detection prospects of single gravitational-wave sources, as well as the detection penalty incurred by employing 
a circular waveform template to search for eccentric signals, and conclude by identifying important avenues 
for future study. 

Subject headings: Gravitational waves - Methods: data analysis - Pulsars: general - 


1. INTRODUCTION 

The observation of extremely compact objects — black 
holes (BHs), neutron stars (NSs), and white dwarfs — and 
the development of a thorough theoretical understanding of 
their nature has been one of the triumphs of modern astro¬ 
physics (Chandrasekhar 1983; Misner et al. 1973; Thome 
1987), but there is still much that we do not understand about 
these exotic objects. The combination of electromagnetic ob¬ 
servations with future detections of gravitational-wave (GW) 
signals will provide key insights into the nature of compact 
objects and the role they play in some of the most energetic 
events in the Universe: gamma-ray bursts, active galactic nu¬ 
clei, quasars, etc. (Hughes & Blandford 2003; Hughes 2009; 
Gebhardt et al. 2000; Soltan 1982; Peterson et al. 2004; Ko- 
rmendy & Richstone 1995; Magorrian et al. 1998; Berger 
2013; Berger et al. 2013; Janka et al. 1999; Lee & Ramirez- 
Ruiz 2007; Metzger & Berger 2012; Piran et al. 2013; Tanvir 
et al. 2013). Several large-scale collaborations are working 
to inaugurate the new field of GW astronomy by targeting 
a wide variety of potential GW sources. These range from 
the mergers of supermassive black hole binaries (SMBHBs), 
which may be used by pulsar timing arrays (PTAs) to probe 
the innermost regions of merging galaxies, to the coalescence 
of NS binaries and stellar mass BHs, which encode important 
information about stellar evolution, galactic nuclei and glob- 
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ular clusters, and are the principle targets for ground-based 
GW detectors. 

In this article, we will focus on a particular type of source 
that is being targeted by PTAs (Poster & Backer 1990). PTAs 
aim to observe GWs in the nanohertz frequency band via the 
accurate timing of millisecond pulsars. There are three major 
PTA collaborations — the European PTA, (EPTA, Kramer & 
Champion 2013), the North American Nanohertz Observatory 
for Gravitational-waves (NANOGrav, McLaughlin 2013) and 
the Parkes PTA (PPTA, Hobbs 2013) in Australia. These three 
collaborations also aim to cooperate as the International PTA 
(IPTA, Manchester & IPTA 2013). 

The sources of interest in this work are individual SMBHBs 
during their early inspiral evolution (Rajagopal & Romani 
1995; Jaffe & Backer 2003; Wyithe & Loeb 2003; Sesana & 
Vecchio 2010; Sesana et al. 2009). Given the nature of these 
systems, i.e., large orbital separations and small local velocity 
of the binary components, we can take the compact objects as 
point-particles without internal dynamics and model the or¬ 
bital evolution of the system using a post-Newtonian expan¬ 
sion (Peters & Mathews 1963; Barack & Cutler 2004; Sesana 
& Vecchio 2010). Purthermore, these events will be observed 
at large orbital separations, where the orbital evolution may be 
more strongly influenced by dynamical interactions with the 
astrophysical environment rather than GW emission. Hence, 
the circularizing influence of the latter may be lessened, al¬ 
lowing for quite large orbital eccentricities at the time of de¬ 
tection. 

There are several mechanisms that could drive the eccen¬ 
tricity evolution of a SMBHB. Por instance, at sub-parsec 
scales a binary formed by a galactic merger may be embed¬ 
ded in a dense stellar environment. As discussed in Sesana 
et al. (2008), if one assumes an isotropic stellar distribu¬ 
tion, the interaction of a star and a SMBHB with semi major 
axis a can have two possible outcomes. Denoting the semi- 
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major axis of the binary formed by the star and the SMBHB 
by encounters with stars with ^ a tend to circularize 
the orbit, whereas those with stars with ^ a tend to in¬ 
crease the eccentricity of the binary. In non-isotropic envi¬ 
ronments, co-rotation of the stellar distribution tends to cir¬ 
cularize the binary. Counter-rotating stars tend to extract an¬ 
gular momentum from the SMBHB, causing the eccentricity 
to grow (Sesana et al. 2011). Several issues still remain to be 
explored regarding the evolution of SMBHBs at sub-parsec 
scales in dense stellar environments, but most models seem to 
favor a growth in orbital eccentricities before these systems 
enter the frequency band of PTAs (Sesana 2010; Roedig & 
Sesana 2012). 

Aside from interactions with stars, the dynamical evolution 
of a SMBHB at sub-parsec orbital separations can also be in¬ 
fluenced by the redistribution of energy and angular momen¬ 
tum between the binary and a self-gravitating disc. Consider 
a gaseous disc co-rotating with a binary, and deflne X = Rt/a, 
where Rt is the distance of the strongest torque on the binary 
as measured from the center of mass, and a is the semi-major 
axis of the binary. Detailed numerical simulations suggest that 
the evolution of the orbital eccentricity of a SMBHB embed¬ 
ded in a circumbinary disc is independent of the mass-ratio of 
the system, but depends sensitively on the location of the inner 
rim of the disc. A, with respect to the binary’s center of mass. 
For 2 < A < 2.5, it is expected that binaries will converge to 
a critical eccentricity value 0.55 < < 0.79. Binaries with 

initial eccentricities e > will undergo a steady decrease in 
eccentricity, whereas binaries with ^ will experience the 
opposite behavior. The larger the separation between the rim 
of the disc and the center of mass of the binary, the longer the 
system will take to attain (Roedig et al. 2011). 

Taking into account these considerations, and the fact that 
uncertainties about the environments of binaries in realistic 
galaxy mergers make binary eccentricity a legitimate possi¬ 
bility, we recently introduced a theoretical framework to ex¬ 
plore in detail the effect of eccentricity for source detection 
of potential PTA sources (Huerta et al. 2015). We now ex¬ 
tend that analysis by introducing novel, accurate and efficient 
pipelines that shed light on the accuracy with which the astro- 
physical parameters of individually resolved eccentric SMB- 
HBs can be reconstructed. This analysis explores the impact 
of eccentricity both in terms of source detection and param¬ 
eter estimation, and presents new statistics to facilitate the 
analysis. Our approach builds on previous Bayesian (Ellis 
2013; Taylor et al. 2014) and frequentist (Babak & Sesana 
2012; Ellis et al. 2012) statistics which have assumed circu¬ 
lar gravitational waveform models, and unlike recent studies 
(Zhu et al. 2015), can recover all binary characteristics in ad¬ 
dition to providing detection statistics. The latter study de- 
flned a frequentist statistic in terms of a harmonic sum over 
the lowest two harmonics, whereas we proceed from the full 
GW strain model of an eccentric binary, producing analytic 
signal models for Bayesian PTA single-source GW searches, 
and a well-motivated frequentist statistic which fully general¬ 
izes that of Babak & Sesana (2012) and Ellis et al. (2012). 

This article is laid out as follows. In Section 2 we briefly 
review the orbital trajectories of eccentric binary systems, and 
how we can analytically solve for the orbital phase at a given 
time. This is followed in Sec. 3 by a description of the ec¬ 
centric gravitational waveforms we use, and in Sec. 4 by our 
model of the perturbations these GWs induce in the times of 
arrival of radio signals from pulsars. The details of our analy¬ 
sis are provided in Sec. 5, followed by the results of Bayesian 



Figure 1. A diagram illustrating the relationship between the various angu¬ 
lar elements in a binary system with orbital eccentricity e, reduced r nass /U , 
and total mass M. The semi-major and semi-minor axes are a and aV 1-e^, 
respectively. If we measure the angles from the moment of periapsis, then 
is the true anomaly, I is the mean anomaly, and u is the eccentric anomaly. 
The auxiliary circle has a radius equal to the orbital semi-major axis. 

and frequentist signal recoveries from simulated datasets in 
Sec. 6. In Sec. 7, we discuss the likely impact of several as¬ 
sumptions that we have made which should be explored fur¬ 
ther in future studies. We flnish with concluding remarks in 
Sec. 8. In the following we adopt units such that G = c = 1. 


2. ECCENTRIC BINARY ORBITS 

We briefly review the Kepler problem and present the gen¬ 
eral approach to analytically solve for the orbit of an eccen¬ 
tric binary, reiterating some of the notation and formalism 
of Yunes et al. (2009), and referring the reader to Goldstein 
(1950) for a more complete discussion. 

We consider a binary system with component 
masses mi and m 2 , total mass M, and a reduced mass 
/i = mim 2 /(mi+m 2 ). The separation vector joining the 
components is deflned in terms of the component position 
vectors by r = ri -r 2 , such that n ^m^f/M and = -mif/M. 
Using (r = |fl,4>) to denote plane polar coordinates for the 
position of one member of the binary with respect to the 
other, the Newtonian Keplerian orbital trajectories of two 
point particles in an eccentric binary system are described by 


r = ( 2(1 -^cosm), 


uj{t — tQ) = l = u-esinu^ 


4>-^o = V = 2arctan 


fl+e 


1/2 


u 

tan - 
2 


( 1 ) 

( 2 ) 

( 3 ) 


where a is the semi-major axis of the orbit, and 0 < ^ < 1 is 
the eccentricity (of a bound orbit). The eccentric anomaly, u, 
is an auxiliary variable with which to parametrize the radial 
and phase coordinates. Given the average angular frequency 
(or mean motion', b: - In jT, where T is the orbital period) 
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and eccentricity of the orbit, we can solve the transcenden¬ 
tal Eq. (2) for at a given time t, where I -luit is 
denoted as the mean anomaly. The eccentric anomaly can 
then be plugged into Eqs. (1) and (3) to give the separation 
and orbital phase (or true anomaly', ^ - ^o) at any point along 
the orbital trajectory. If we assume that time and phase are 
measured from the moment of periapsis, then the constants of 
integration and T>o can be set to zero. All of these angular 
quantities are shown diagrammatically for an example orbital 
ellipse in Eig. 1. 

The flux of energy and angular-momentum carried away 
from the system by GWs depend on the eccentricity and the 
Keplerian mean orbital frequency, F. Once the binary evo¬ 
lution is driven solely by GW emission, these co-evolve as 
(Peters 1964) 


where 


F{e) 


(4) 

Fieo) 

V <7(e) J ’ 

^12/19 r 

■ 121 


cr(e)= 2 

\-e^ 


(5) 


and eo is deflned as the eccentricity of the system at some 
earlier reference epoch of the binary evolution. 

The frequency F can be regarded as the instantaneous mean 
orbital frequency. Eor GW-dominated orbital evolution, it co¬ 
evolves with the eccentricity according to the coupled differ¬ 
ential equations (Peters 1964) 


&F 

dt 

de 

dt 


48 

57rA4^ 

304 

~15M 


(27rMFy^^ 


1 + 73.2+37 4 
(l_^2)7/2 ’ 

1 + ^e^ 

^ ^ 304 ^ 


( 1 -^ 2 ) 5 / 2 ’ 


( 6 ) 


where A4 = {mimfy"^^l(m\ +m 2 )^/^ is the binary chirp mass. 

Gravitational waveform templates describing the emission 
from inspiraling binary systems depend on trigonometric 
functions of the orbital phase. Eor circular systems the re¬ 
lationship between orbital frequency, time, and phase is sim¬ 
ple: we have ^ = 2tt j Fit)dt, where Fit) is the Keplerian 
orbital-frequency (half of the dominant quadrupole GW fre¬ 
quency) which evolves according to Eq. (6) with e = 0. How¬ 
ever the situation is rather more complicated for eccentric sys¬ 
tems. The phase is related via an arctangent to the eccentric 
anomaly, which is then related to the mean anomaly (and thus 
the mean angular frequency uj = IttF) via a transcendental 
equation. The so-called Kepler problem refers to the histori¬ 
cal difficulty in finding solutions to the transcendental equa¬ 
tion in Eq. (2) and thus being able to express the orbital phase 
in terms of the mean anomaly. We do so using the well known 
Eourier analysis of the Kepler problem. Eor full details of the 
calculation see Watson (1995). Using elementary properties 
of elliptic curves and Bessel functions, the results are 


2 

cos$ = -^+-(l-^^)y^ Jnine)co^inl), (7) 

e 

n=\ 

oo 

sin$ = (l-e^)*/^y^[/„_i(ne)-y„+i(ne)]sin(n/). (8) 

n=l 



Figure 2. The minimum number of harmonics required for the Fourier so¬ 
lution of cos as a function of I (mean anomaly) to maintain accuracy with 
the numerical solution. We demand that the overlap of the Fourier solution 
and numerical solution, as determined by the normalized scalar product of 
the two solution vectors, is > 99.999% over of mean anomaly. 

Spiraling binary systems in terms of the mean orbital fre¬ 
quency. Equations (7) and (8) can be immediately used to 
construct these waveforms. However, by setting a required 
tolerance on the accuracy of sin 4> and cos 4> for a given eccen¬ 
tricity, we can truncate the infinite summations (Pierro et al. 
2001; Yunes et al. 2009) to accelerate calculations. We inves¬ 
tigate the minimum number of terms required for the Eourier 
series expansion of cos in Eq. (7) to maintain accuracy with 
the exact numerical solution of Eqs. (l)-(3), by demanding 
that the error in the two solutions (determined by the normal¬ 
ized scalar product between the two solution vectors) is less 
than 0.001% over 27r of mean-anomaly. The results are shown 
in Eig. 2, where we see that ^100 terms in the summation 
are necessary to maintain accuracy up to ^ = 0.9, however the 
required number of terms dramatically increases beyond 0.9, 
exceeding 10^ at ^ = 0.99. 

Although systems with high residual eccentricity (> 0.9) 
in the sub-parsec GW inspiral regime may exist, they are 
by no means expected to be common. Unequal mass sys¬ 
tems with q < 0.25 may retain ^ > 0.9 into the PTA band 
(Khan et al. 2012), but we are unlikely to detect their weaker 
gravitational-wave emission, so we focus here on the more 
probable case of a detectable signal from a comparable mass 
binary. Comparable mass binaries in isolated galaxy sim¬ 
ulations exhibit e < 0.95 when they transition from stel¬ 
lar hardening to gravitational-wave-dominated evolution, al¬ 
though preliminary merger simulations can produce binaries 
with larger eccentricity (Vasiliev et al. 2015). However, since 
gravitational-wave emission is well known to decrease ec¬ 
centricity (Peters & Mathews 1963), we believe the assump¬ 
tion that most detectable systems will likely have ^ < 0.9 in 
the PTA band is astrophysically well-motivated, in addition 
to simplifying things computationally. Ultimately, the range 
of orbital separations at which the transition between stellar 
hardening and radiation-reaction occurs in real galaxies is a 
matter of debate (along with the range of possible binary ec¬ 
centricities in the PTA band), and may only be resolved with 
pulsar-timing measurements. Hence, in the following we re¬ 
strict our attention to systems with eccentricity below 0.9, in 
which regime highly accurate waveforms require the inclu¬ 
sion of fewer than 100 Eourier terms. 


With these trigonometric functions of the orbital phase, we 
can now construct gravitational waveforms for eccentric in- 


3. ECCENTRIC TIME-DOMAIN WAVEEORMS 
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In the transverse-traceless gauge the GW-tensor can be 
written as a linear superposition of “plus” and “cross” po¬ 
larization modes, with 
/z|+x}, and basis-tensors, 

habit A) = + ( 9 ) 

where ft is defined as the direction of GW propagation. 

We employ the Peters-Mathews waveforms (Peters & 
Mathews 1963) given by Barack & Cutler (2004), which make 
use of the Fourier analysis of the Kepler problem to give the 
following analytic expressions for h+ and : 

h+(t) = -(1 + cos^ i){an cos( 27 ) - bn sin( 27 )] 

n 

-h(l-COS^ L)Cn^ 

hx{t) = 2 COS L[bn cos( 27) -i- Un sin(27)] , (10) 


associated polarization-amplitudes, 
'’^'^’^^(^2), such that 


'=^ab 


where 

an [/„_2(^^)-2^/„_i(^^)-r(2/^)/„(^^) 

-t- 2eJn+\ {ne)-Jn+ 2 {ne)^ cos[^/(0], 

bn =-nCuj'^^^^/l^[Jn- 2 {ne)- 2 Jn{ne)+Jn+ 2 {ne)] sin[^/(0], 
= 2Cuj^^^Jn(ne)cos[nl(t)]. ( 11 ) 

The amplitude parameter is defined as C = jD^, where 
is the luminosity distance of the binary, and uj = 2 tvF. The 
mean anomaly is lit) = /o+27r Fit')dt' (where /q is the mean 
anomaly at ^o); 7 is an azimutlial angle nieasuring the direc- 
tion of pericenter with respect to v = (^2+Lcos ^)/Vl-cos^^; 
and L is the binary orbital inclination angle, defined by 
cos i = -L'(l. In the following, F and M. refer to the observed 
redshifted values, such that /y = F(1 -Hz) and Afr = Af/(l+z), 
where /v and Mr are rest frame values, and z is the cosmo¬ 
logical redshift of the binary. 

An important feature to emphasize here is that eccentric bi¬ 
naries do not radiate monochromatic GWs, but rather emit 
a spectrum of frequencies which are harmonics of the mean 
orbital frequency. Given that /o(0) = 1 and /„>o(0) = 0, it is 
immediately obvious from Eqs. (10) and (11) that e = Q wave¬ 
forms will only include the ^ = 2 harmonic of the binary’s 
mean orbital frequency. This is the usual result that the GW 
frequency of emission from circular binaries is twice the or¬ 
bital frequency. 

To construct the polarization basis tensors, we define a 


right-handed basis triad in terms of where h = -tl, 

p^(nxL)l\n X L\ and q = pxh. The vectors comprising the 
basis triad are explicitly 

=(sin6>cos0, sin6>sin0,cos6>), (12) 

p = (cos cos 0 cos 0 - sin sin 0, 

cos t/; cos 6> sin 0 -h sin ip cos 0, - cos ^psinO), (13) 

q = (sin pj cos 0 cos (j) + cos pj sin 0, 

sin pj cos Osincp- cos cos 0, - sin pj sin 6>), (14) 


where iO^cp) = (7r/2-DEC,RA) denotes the sky-location of 
the binary in spherical polar coordinates, and pj corresponds 
to the angle between p and the line of constant azimuth when 
the orbit is viewed from the origin of our coordinate system. 
These angles are shown diagrammatically in Fig. 3. The vec¬ 
tors p and q lie in the plane that is transverse to the direction 


L 



Figure 3. A diagram illustrating the geometry of an eccentric SMBHB with 
respect to the angles of our coordinate system. The unit vector pointing to the 
binary is n = -Cl, with spherical-polar coordinates {6 = 7r/2 -DEC, 4> = RA}. 
The binary orbital inclination angle is defined by cos i = L-n, where L is a 
unit vector pointing along the binary’s orbital angular momentum. The GW 
polarization basis tensors are defined in the plane transverse to the direction 
of propagation, in terms of the unit vectors p = (hxL)/\hxL\ and q = pxn, 
where {n,p,q} define a right-handed basis triad. The vector p lies along the 
major axis of the projected ellipse as seen from the origin of the coordinate 
system. The GW polarization angle 'ip is defined as the angle between p and 
the line of constant azimuth. This diagram is a modified version of Fig. 1 in 
Apostolatos et al. (1994). 

of GW propagation, and are used to construct basis tensors as 
follows: 


e'lb=PaPb-qaqb, ( 15 ) 

^ab = P‘>qb + qaPb- ( 16 ) 

4. PULSAR TIMING RESIDUALS INDUCED BY AN 
ECCENTRIC BINARY 

As a GW transits across the line of sight between a pulsar 
and the Earth, it creates a perturbation in the space-time met¬ 
ric which causes a change in the proper separation between 
the Earth and the pulsar. This in turn leads to a shift in the per¬ 
ceived pulsar rotational frequency. The fractional frequency 
shift of a signal from a pulsar in the direction of unit vector 
u, induced by the passage of a single GW propagating in the 
direction of Q is (Anholm et al. 2009; Book & Elanagan 2011) 

z(tP)=lP^-^—Ahab(tP), ( 17 ) 

2 1+D-m 

where Ahab = habiteM)-habitpM) is the difference in the 
metric perturbation evaluated at time 4 when the GW passed 
the solar system barycenter (SSB) and time tp when the GW 
passed the pulsar. Froni simple geometrical arguments, we 
can write = 4 - L(l + (2 • fi), where L is the distance to the 
pulsar. The integrated effect of this GW-induced redshift over 
the total observing time of the pulsar leads to an offset be¬ 
tween the expected and the observed pulse TOA: 

s{t)= [ zit')dt'. (18) 

Jo 

The expected pulse TOA is computed from a deterministic 
timing model which characterizes a pulsar’s astrometric and 
spin properties. This model is refined over many observations 
to give an accurate prediction of the pulse arrival times. The 
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difference between the measured TOAs and those predicted 
by the best-fit deterministic timing-model are the timing resid¬ 
uals. In addition to any GW signals, these residuals encode 
the infiuence of noise processes and all unmodelled phenom¬ 
ena which affect pulsar TOAs. The pulsar timing residuals 
induced by a single GW source can be written as 

sit, h) = F^ih)/\s^it) + F^ih)/\s,,it), (19) 

where A = {-h, x|,^ A^a(0 = SAite) - SAitp), with ^a( 0 = 
/o hAit')dt', and F^iFt) are antenna pattern response functions 
encoding the geometrical sensitivity of a particular pulsar to a 
propagating GW, defined as 

( 20 ) 

and corresponding to the contraction of the pulsar-timing im¬ 
pulse response function with the GW polarization basis ten¬ 
sors. 

The form of ^-aCO can be computed analytically by as¬ 
suming that the binary’s mean orbital frequency and eccen¬ 
tricity remain constant over the total timespan of our obser¬ 
vations of a given pulsar. More specifically, we must as¬ 
sume no binary evolution over the Earth term timing baseline, 
[tefe + T], and also the corresponding timing baseline of the 
pulsar term, {te-Li\-\-^-u),te + T-Lil + ^-u)], where T is 
0(10 years). ^ Therefore, time only appears in the definition 
of the mean anomaly as a linear parameter, such that lit) = 
Iq + 2ti f^^Fit')&t' = iQ + luFit -tf), which allows cos[^/(0] 
and sin[^/(0] in Eq. (11) to be trivially integrated to give the 
plus/cross residuals: 

^+(0 = E -(1 -h cos^ i){cin cos( 27 ) - Sn sin( 27 )] 

n 

-l-(l-COS^ L)Cn^ 

5'x(0 =^2cos^[^„cos(27)-l-i7„sin(27)], ( 21 ) 

n 

where 

an= [Jn-2{ne)-2eJn-\{ne) + {2/ri)Jn{ne) 

+ 2 eJn+\{ne)-Jn+ 2 {ne)] sin[n/(0] 

= v ^sin[ ^/( 0 ], 

Sn = -e^ Un -2 ine)-2Jn ine)+Jn +2 (/t^)] cos {nl (01 

= X4COS[^/(0], 

Cn = i2/n)(luj~^^^Jnine)^m{nlit)] 

= v,^sin[^/(0], (22) 

and the quantities {xa^,xs^,Xcf\ are defined for later conve¬ 
nience. 

We can now analyze the harmonic content of the vari¬ 
ance of the residuals from both plus and cross polarizations, 
which is computed over one period of binary elliptical mo¬ 
tion (/ = {0,27r}) and over cos ^, 7 . Clearly averaging over a 
single (or any non-zero integer) period of orbital motion is 
only an approximation, since our pulsar-timing observations 
are highly unlikely to span an integer number of orbital peri¬ 
ods or GW cycles. Nevertheless we carry out this calculation 
since it illuminates certain features of the harmonic content 

^ The binary’s mean orbital frequency and eccentricity do evolve non- 
negligibly over the light travel time between the Earth and the pulsar, 
G(1000 years)-G( 10000 years). This effect is easily included in our signal 
model, however in the rest of this paper we consider only the Earth term. 


of the GW signal from eccentric SMBHBs. We employ the 
following relations when averaging over the mean anomaly: 

pin 

/ dl sininl) cosin'1) = 0, y n,n', (23) 

Jo 

[ dl sminl)smin'l)= I (24) 

Jo [tt, if n = n, 

where n,n' > 1 , and the last equation is also true for cosine 
functions. Given that the induced residuals are zero-mean 
over integers of the binary orbital period, the resulting vari¬ 
ance of the residuals is 

= (25) 

n 

where 

(4)« = ^(7+7)+^^t 

(4)« = ^bl+4)- (26) 

The value of {s\)n for several binary eccentricities is shown 
in the left panel of Eig. 4. At each eccentricity, the contribu¬ 
tion of each harmonic to the variance of the residuals is nor¬ 
malized with respect to the largest contribution. In the right 
panel of Eig. 4 we show the fraction of the total variance of 
the plus-component timing residuals contributed by the dom¬ 
inant harmonic, which switches from ^ = 2 in the 0 < ^ < 0.4 
range to ^ = 1 beyond e ^ 0.4. 

Eor the remainder of this paper we will present results from 
investigations with the Earth term of the GW-induced timing 
residuals. The signal model in Eq. (21) is general, and can 
be used to compute both Earth and pulsar terms, modulo the 
assumption of binary non-evolution over typical pulsar tim¬ 
ing baselines. However, including the pulsar term requires 
either precise knowledge of the individual pulsar distances, 
or the distances to be searched or marginalized over (Ellis 
2013; Taylor et al. 2014). This search over distance brings 
its own challenges since the likelihood is highly sensitive to 
small changes in the sampled distance around the true value, 
and can lead to inefficient sampling. We defer considerations 
of the pulsar term to future work, but will briefiy consider its 
infiuence in Sec. 7. Eurthermore, for the most extreme combi¬ 
nations of binary mass, eccentricity, and orbital frequency, the 
system may exhibit frequency chirping and orbital circulariza¬ 
tion during typical pulsar-timing observation timespans, ren¬ 
dering the assumption of non-evolution invalid. We explore 
these issues in Sec. 7 amid suggestions for future directions. 

Related to these two issues are the fact that in general we 
would also need to consider evolution of the direction of peri- 
center, 7, and orbital plane precession from spin-orbit cou¬ 
pling. Evolution of the direction of pericenter can occur even 
for circular binary systems composed of non-spinning black 
holes, leading to phase shifts and recovery bias in the orbital 
frequency if it is not considered. However, as discussed in 
Sesana & Vecchio (2010), these factors can be safely ignored 
over typical PTA observation timespans. In Eig. 5 we show 
exclusion regions in {M = (mi +m 2 ),F,e} parameter space, 
where pericenter direction evolution leads to a bias in the 
orbital frequency which is greater than the typical PTA fre¬ 
quency resolution of l/T for a 10 year observation timespan 
(Sesana & Vecchio 2010). The excluded regions correspond 
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Figure 4. (Left): The contribution of each harmonic of the orbital frequency to the variance of the plus-component timing residuals. At each eccentricity 
we normalize the contributions from each harmonic with respect to the maximum contribution. The only contribution for circular binaries is from the second 
harmonic (black star and line, slightly offset from n = 2 for ease of viewing). At higher eccentricities (e = 0.5,0.9) the contribution is spread into a spectrum 
of higher harmonics, but is dominated by the fundamental harmonic. (Right): The fraction of the total variance contributed by the dominant harmonic, h, as a 
function of eccentricity. As in the left panel, n labels the harmonic of the binary mean orbital frequency. In the range 0 < ^ < 0.4 the second harmonic dominates, 
whilst beyond e ~ 0.4 the fundamental harmonic dominates the variance of the induced timing residuals. 
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Figure 5. Exclusion regions in binary eccentricity and orbital frequency as 
a function of binary total mass, corresponding to parameter combinations 
where unmodelled evolution of binary pericenter direction causes a bias in 
orbital frequency recovery which could be resolved by 10 years of PTA ob¬ 
servations, A/ = 1/r = 3.2 nHz. 


to systems with very high total mass and eccentricity, and or¬ 
bital frequencies beyond the region of peak PTA sensitivity. 
Hence, we ignore this effect here and consider only {F,e} 
evolution in Sec. 7, but information from these additional ef¬ 
fects may allow the individual binary component masses, and 
possibly their spin, to be constrained (Mingarelli et al. 2012). 
Additionally, these effects are likely to be highly important 
when tracing the binary evolution back by thousands of years 
to the pulsar term. 

5. SIMULATED DATASETS AND ANALYSIS 

Eor our proof-of-principle study of an eccentric single¬ 
source pipeline, we consider two types of PTA datasets. In 
our Type I array, we consider the 36 pulsars from the IPTA 
mock data challenge.^ They are timed to 100 ns precision 


over a timing baseline of 10 years, with observations carried 
out every 4 weeks. This array is obviously idealized, how¬ 
ever the generalization to more realistic observing schedules 
and pulsar noise properties does not require modifications to 
our pipeline since it is constructed in the time-domain, and is 
shielded from Eourier domain spectral leakage caused by red 
timing noise or irregular sampling. The Bayesian pipeline can 
be trivially incorporated into a more general pipeline which 
simultaneously estimates pulsar noise properties and other 
stochastic signals. The Type I datasets will serve as the ideal 
observing scenario to test for any systematic errors in our sig¬ 
nal construction which are separate from observing practical¬ 
ities, and will also be used for brief analyses of the influence 
of binary eccentricity on circular- or eccentric-model signal- 
to-noise ratios (SNRs). 

To emulate more realistic observing schedules and pulsar 
noise properties, we also construct Type II datasets using the 
actual epochs of observation and noise properties of the 18 
pulsars that were used by the NANOGrav collaboration to 
place astrophysical constraints on the nanohertz GW back¬ 
ground (The NANOGrav Collaboration et al. 2015; Arzou¬ 
manian et al. 2015). These pulsars suffer from irregular sam¬ 
pling, different timing baselines (the longest is ^ 9 years), 
heteroscedastic TOA measurement errors, and, in some cases, 
intrinsic pulsar spin noise. These Type II arrays will be used 
for our Bayesian studies of the penalties arising from assum¬ 
ing a circular binary model when analyzing data having an 
eccentric signal, and also when estimating the precision with 
which current PTAs can estimate binary parameters. 

We use the simulation routines within libstempo,^ a python 
wrapper for the pulsar-timing software package TEMP02 
(Hobbs et al. 2006; Edwards et al. 2006). Eor a fiducial 
source, we are only interested in sensible binary parameters 
which will illustrate the efficacy of the search pipeline. We 
follow Ellis (2013); Taylor et al. (2014) by considering a 
source with the following characteristics: = = 

5 nHz, (j) = 0.95,= 2.17, ^ = 1.57, /q = 0.99, = 1.26 ,7 = 0.5}, 
and a luminosity distance scaled to meet a required optimal 


^ http://www.ipta4gw.org/?page_id=89 


^ http://vallis.github.io/libstempo/ 
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SNR. The definitions of optimal and matched-filtering SNR 
follow from Finn (2001). 

The binary parameter space is searched using a python 
wrapper (Buchner et al. 2014) to the nested sampling pack¬ 
age MultiNest (Feroz et al. 2009; Feroz & Hobson 2008; 
Feroz et al. 2013), and we have cross-checked our results with 
a sampler utilizing advanced Markov chain Monte Carlo tech¬ 
niques.^^ The product of these analyses are samples from 
the posterior probability distribution of the signal parame¬ 
ters space, allowing us to quantify the measurement preci¬ 
sion of parameters based on the Bayesian credible regions, 
and also permitting model selection via computation of com¬ 
peting models’ Bayesian evidence. The priors for the sig¬ 
nal parameters are as follows: log^Q(7W/M0) G f/[7,10], 
logio(D^/Mpc) G f/[0,4], logio(F/Hz) G f/[-9.3 -6.0], ^ G 
f/[0,0.9], (j) G f/[0,27r], cos6> G f/[-l,l], cos^ G f/[-l,l], 

G f/[0,7r], 7 G t/[0,7r], /q G f/[0,27r]. 

Full details of Bayesian inference in the context of PTAs 
can be found in van Haasteren et al. (2009); van Haasteren & 
Vallisneri (2014); Arzoumanian et al. (2015), and for details 
of how Bayesian searches for continuous GWs are carried out 
see Ellis (2013); Taylor et al. (2014). 


5.1. Eccentric statistic 

Equation (21) provides the appropriate signal model to use 
when we wish to map out the posterior distribution of the 
entire signal parameter space, and also if we were to simul¬ 
taneously search for continuous GW sources in addition to 
stochastic signal or noise processes. However, we can also 
construct a fixed-noise frequentist statistic for eccentric bi¬ 
nary systems. 

We wish to construct a form of the statistic (Babak & 
Sesana 2012; Ellis et al. 2012) which can be applied to GW 
signals from binaries with arbitrary eccentricity. In practice, 
as in the rest of this paper, we only consider systems with 
e G [0,0.9]. The statistic as it is constructed in Ellis et al. 
(2012) is a maximum-likelihood estimator of the source’s 
sky-location and orbital frequency, and requires that the ex¬ 
pression for the induced residuals be rearranged into a form 
which permits maximization of the likelihood-ratio over the 
coefficients of a set of time-dependent basis-functions. The 
likelihood-ratio. A, is defined as the ratio of the likelihood of 
the data in a model which includes a signal to the noise-only 
null hypothesis: 


lnA=ln 


>C(s|(5t) 

C(0\6i) 


=((5t|s)--(s|s). 


(27) 


We extend the statistic by rewriting the Earth term resid¬ 
uals (in a single pulsar) given by Eqs. (21) and (22) as: 


s(i) = Y,nW‘, 


(28) 


i=l 


https://github.com/jellislS/PTMCMCSampler 


where, 

ivi =( [-(l-fcos^^) cos(27)cos(2'0)-h2cos^sin(27)sin(2'0)] , 
ZC 2 =C [(l+cos^^) sin(27)cos(2'0)-t-2cos^cos(27)sin(2'0)] , 
zi^ =( [(1-cos^ ^)cos(2'0)] , 

=C [ (1 + cos^ l) cos( 27) sin(2'0) -h 2 cos l sm(2^) cos(2'0)] , 
zcs =( [-(l-tcos^ sin(27)sin(2'0)-l-2cos^cos(27)cos(2'0)] , 
^6 =C [-(1 -cos^ l) sin(2'0)] , (29) 


sm[nuj(t-to)-\-nlo]^ 

n 

co^[nuj{t-t{)) + nl{)\, 

n 

-^41 sin[^C(;(r-ro) + ^/o]5 

n 

^m\nuj{t- tQ) + nlQ\, 

n 

yyS _px ^ co^[nuj{t - to )+ hIq] , 

n 

y^6 _ sin[nuj(t-to)-\-nlo], (30) 

n 

and we adapt the number of terms in these summations based 
on the binary eccentricity. This is the same adaptation as dis¬ 
cussed in the previous section for the Bayesian analysis. 

The antenna pattern functions F^(Sl) are related to F^ijl) 

by, 

(F^\_( cos( 2V^) -sin(2V^)\ (F^\ ...x 

) \^sin(27/;) cos(27/;)y \F^ j ' ^ 

The coefficients zvt are a function of extrinsic source parame¬ 
ters {C, L, 7}, whilst the time-dependent basis-functions W' 
are a function of intrinsic source parameters {F, 0^(t),e^lo}. 
Hence, the full PTA signal template can be written as: 


6 

s(f) = ^«^,WV), (32) 

i=l 


where, 


'Wj(r)' 

>V'(r) 


(33) 


and yVj(t) denotes the quantity W' defined by Eqs (30) for 
pulsar j. Inserting Eq. (32) into Eq. (27) and using Einstein 
summation convention, we have 

In A = zvi zvj , (34) 

where W = ((5t|W0 and = (W^|W^). By maximizing the 
log-likelihood ratio over the amplitude coefficients, zvi, we get 
their maximum-likelihood values: 


ici = MijN\ (35) 

where Mfj = Substituting these coefficients back 


Through practical experience we find that the inverted matrix has greater 
numerical stability at low eccentricity (e < 0.05) when a Moore-Penrose 
pseudoinverse is used, with a typical singular value cutoff of ~ 10“^^. 
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into the expression for In A gives the eccentric Te statistic: 

(36) 

The procedure to estimate the maximum likelihood values 
of all of the signal parameters is as follows: 

• We find the local maxima of the Te statistic in the space 
of intrinsic parameters via a straightforward function 
maximization, or we can map out the posterior distri¬ 
bution of the semi-maximized parameter space with a 
stochastic sampler, and determine the maximum likeli¬ 
hood point from the resulting chain. 

• The intrinsic parameters which maximize the Te statis¬ 
tic can be used to compute the quantities Mij and 
which are combined to determine the maximum likeli¬ 
hood coefficients, zvi, via Eq. (35). 

• From these coefficients, we obtain a maximum likeli¬ 
hood estimate of the physical extrinsic parameters, as 
described below. 

There are six wi parameters, but these are functions of only 
four physical extrinsic parameters and so not all combinations 
of correspond to physical systems. However, we can ob¬ 
tain extrinsic parameter estimates from estimates of the 
following (Cornish & Porter 2007). We define 


from a single pulsar in our array. The results are shown for an 
e = 0 and ^ = 0.5 binary signal in Fig. 6, where we see that the 
region of credible residuals (enclosed within red dashed lines) 
tracks the main features in the raw post-fit residuals, and cor¬ 
rectly interprets high frequency behavior around MJD 55100 
in the right panel (e = 0.5) as binary periapsis. In Fig. 6 we 
also show the deviation of the recovered residuals from the 
true injected residuals, where the envelope of credible residu¬ 
als encompasses the line of zero offset. This shows that, even 
in this high SNR case, any systematic bias from the adaptation 
of the number of harmonics is very small, and our Bayesian 
pipeline is robustly recovering the signal characteristics. 

We also use our samplers to map out the Te statistic dis¬ 
tribution over the intrinsic parameter space. From the chain 
of sampled points we determine the maximum-likelihood in¬ 
trinsic parameters, which are then used to construct via 
Eq. (35). Having the maximized zi^i and corresponding 
we now compute the maximum-likelihood timing residuals 
induced by the GWs from an eccentric binary. The results for 
the e = 0 and ^ = 0.5 binary signals are shown in Fig. 6, where 
the maximum-likelihood GW-induced post-fit residuals are 
overlaid as black dashed lines on top of the raw post-fit resid¬ 
uals from pulsar J0030-I-0451, showing excellent tracking of 
the residual behavior and good agreement with the Bayesian 
recovery. Note that these maximum-likelihood residuals are 
offset by -I-O.l /is for ease of viewing. 


A+=\/ {zvi—zv^^ + , 

A X=v izvi+zv^'^ + {zV2-'^d^ \J (zvi-zvs)'^ + , (37) 

and ! - 

A=A^ + ^Al-A\. (38) 

By employing the quantities {A+,Ax,A} we can map from 
^/g[i, 6] to {C, V^,7} with the following manipulations: 



tan(27/;) = 


A^zvi-A+zvs 
Ax '^4 +A+ZV2 ’ 


tan( 27 ) = 


Ay^zVj-A+zvs 
A+zV/\. + Ax zx^i 


In the following we treat the eccentric Te statistic as 
lihood function and map out the posterior probability 
bution of the semi-maximized signal parameter space. 


(39) 

a like- 
distri- 


6. RESUFTS 
6.1. Efficacy of pipelines 

We accelerate the generation of templates for the GW- 
induced residuals by making the number of waveform har¬ 
monics adapt based on the current proposed eccentricity. As 
discussed in Sec. 2, the number of harmonics to adequately 
describe a binary with ^ = 0.5 is ^ 10, whilst for ^ = 0.9 it is 
^ 100. Adaptation of the number of harmonics avoids tem¬ 
plate generation being the main computational bottleneck in 
our pipeline. 

As a first illustration of the efficacy of our pipelines, we 
inject GW signals with SNR = 20 into noisy Type I datasets, 
and analyze the data with our Bayesian and frequentist statis¬ 
tics. We overlay the 95% envelope of Bayesian credible post- 
fit GW-induced residuals on top of the raw post-fit residuals 


6.2. Detection prospects & parameter precision 

One might expect that distinctive high-frequency features 
due to periapsis passage (such as seen in Fig. 6) may improve 
the prospects for detection. We investigate this by computing 
the optimal SNR for a binary with varying orbital frequency, 
and a PTA timing baseline of 10 years in Type I data. We draw 
the angular waveform parameters randomly and average over 
the resulting SNRs. The result of this procedure as a func¬ 
tion of binary eccentricity is shown in Fig. 7, where we see 
a transition in behavior as the binary orbital frequency moves 
through the most sensitive location in the pulsar-timing band. 
From theoretical calculations and analysis of real data (Moore 
et al. 2015; Yardley et al. 2010; Arzoumanian et al. 2014), we 
expect the region of peak PTA sensitivity to a continuous GW 
to be at a GW frequency of ^ IfT -2 fT. Sensitivity is in¬ 
hibited at lower frequencies by fitting of the pulsar quadratic 
spindown parameters in its timing-model, and higher frequen¬ 
cies are dominated by white TOA measurement errors. For 
^ = 0 binary signals in this simulated PTA, this peak corre¬ 
sponds to an orbital frequency of ^ 1.6-3.2 nHz. In Fig. 7 
we see that at higher eccentricities the SNR is enhanced when 
the injected orbital frequency lies below 1 nHz, and dimin¬ 
ished when it lies above 5 nHz. We can make sense of this 
by recalling the spectral decomposition of the variance of the 
GW-induced residuals shown in Fig. 4, where as the eccentric¬ 
ity is increased the variance is distributed amongst higher har¬ 
monics of the orbital frequency. For systems with F < 1 nHz 
this will enhance the SNR since power in the residual vari¬ 
ance is shifted into the region of peak PTA sensitivity, while 
for systems with F >5 nHz this diminishes the SNR since the 
power in the residual variance is distributed into higher, less 
sensitive frequencies of the PTA band. 

With the approximate scaling behavior of SNR with signal 
eccentricity established, we now investigate how SNR maps 
to the Bayes factor of a signaHnoise model versus a noise 
model alone. Bayesian model selection is actually carried out 
by computing the posterior odds ratio, which is the ratio of 
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Figure 6. The post-fit residuals of pulsar J0030+0451 for simulated Type 1 data are shown in the upper portions of both panels as blue points with associated 
error bars. The left panel corresponds to an injected GW signal from a circular (^ = 0.0) binary, while the right panel corresponds to an injected GW signal from 
an e = 0.5 binary. (Upper): The boundaries of the 95% credible envelope of post-fit residuals induced by the GWs are shown as red dashed lines, while the 
residuals corresponding to the mean signal parameters are shown as solid black. These GW residuals are computed from the parameter posterior PDFs returned 
by Bayesian analysis of the simulated data, and then projected to post-fit values (Demorest et al. 2013). The black dashed line shows the maximum likelihood 
post-fit residuals returned by an eccentric J>-statistic (see Sec. 5.1) analysis (residuals are offset by -l-O.l piS for ease of viewing). (Lower): The offset of the 
reconstructed GW-induced residuals from the injected residuals is shown, where all lines correspond to the same cases as the upper panels. The boundaries of 
the 95% Bayesian credible envelope of post-fit residuals encompasses A = 0, which is a good indicator of the robustness of the pipeline. 



Figure 7. Normalized optimal SNR of a single source as a function of the 
binary eccentricity for a PTA timing baseline of 10 years (Type I data). Only 
the Earth-term component is considered. Each curve corresponds to a differ¬ 
ent choice of binary orbital frequency, and is computed by averaging the SNR 
over all waveform angular parameters. For reference, the GW frequency of 
greatest sensitivity in this Type 1 pulsar array is ~ 5 nHz. 
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Figure 8. Bayes factors for a signal+noise model versus a noise model alone 
in Type II datasets with varying SNR injections. The injected binary param¬ 
eters are the fiducial values given at the start of Sec. 5. The solid blue line 
shows the results for the full eccentric Bayesian pipeline, while the dashed 
green line shows the results for searches over the semi-maximized signal pa¬ 
rameter space (intrinsic parameters) in the eccentric statistic. Red lines in 
the inset figure show the SNR at which each technique reaches a Bayes factor 
of 100. 
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Figure 9. Expected binary parameter measurement precisions from varying SNR injections into Type II datasets. The values plotted are the full width of the 
respective Bayesian credible regions. The hatched regions correspond to the parameter’s prior boundary. The injections are noiseless (but with pulsar noise 
characteristics modeled in the likelihood functions) and carried out at F = 5 nHz, ^ = 0.5. The top row corresponds to precisions deduced from the width of 
Bayesian credible regions produced by our Bayesian pipeline analysis. The bottom row corresponds to precisions deduced from the width of Bayesian credible 
regions produced by mapping the posterior distribution of the intrinsic parameter space with the eccentric statistic. 


competing model evidences (Bayes factor, B) multiplied by 
the prior odds ratio of each model. However, in the following 
we treat the latter quantity as being unity since in real searches 
we can not judge the a priori odds of a signal being in our 
data. Assuming fixed noise properties, the computation of 
the Bayes factor follows by integrating the likelihood ratio in 
Eq. (27) over the signal parameter space. We judge a model 
to be favored over another if the Bayes factor exceeds 100 
(lnS>4.6). 

We inject varying SNR signals into noiseless Type II 
datasets with fiducial parameters and ^ = 0.5. This eccentricity 
is a compromise between being a moderate value in our range 
of exploration, and (as seen later) where the discrimination 
between eccentric versus circular signal models is greatest. 
The injections are noiseless so that we can avoid the need to 
perform a large program of injections to average over noise 
realizations (Nissanke et al. 2010; Cornish 2010). There are 
some reservations over this approach in the low SNR regime 
(Vallisneri 2011), however it is nevertheless correct at high 
SNR, so that we can consider the conclusions drawn here as 
optimistic but indicative of general trends. Importantly, in the 
following we verified that the peak of the recovered posterior 
distributions matched the injected signal parameters, which 
should be the case when the datasets are noiseless. This con¬ 
firms that our techniques do not suffer as a result of the ir¬ 
regular sampling and heteroscedastic uncertainties associated 
with real data. 

Figure 8 shows the growth of Bayes factors favoring a sig- 
nal-fnoise model over a noise model alone in both the full sig¬ 
nal parameter space (Bayesian pipeline) and intrinsic param¬ 
eter space (eccentric statistic). Since the statistic is al¬ 
ready maximized over half of the full signal parameter space, 
it has a lower search dimensionality than the full Bayesian 
pipeline and thus receives less of an Occam penalty. As seen 
in the inset of Fig. 8 the statistic reaches a Bayesian de¬ 
tection threshold at SNR ~ 5 whilst the full Bayesian pipeline 
does so at SNR ^ 7. 


An important question associated with GW detection is 
whether a threshold signal will be associated with any mean¬ 
ingful parameter measurement precisions. We address this 
issue by analyzing the widths of the {68%,95%,99.7%} 
Bayesian credible regions with respect to the prior widths 
for binary orbital frequency, eccentricity, and sky location, at 
varying SNR. Our datasets are again of Type II with fiducial 
parameters and e = 0.5. The results are shown in Fig. 9, with 
measurement precisions obtained using the Bayesian pipeline 
along the top row, and precisions obtained with the eccentric 
IFe Statistic along the bottom row. At high SNR the precisions 
for logioF and e obey a 1/SNR scaling, and the sky location 
(being a compound of two parameters) obeys a 1 /SNR^ scal¬ 
ing. Both techniques perform comparably for SNR > 8, how¬ 
ever if we look at the width of the 99.7% credible region, we 
see that the eccentric statistic begins to update our prior 
knowledge of the parameter space at SNR > 5, whilst this 
happens at SNR > 7 for the full Bayesian pipeline. Hence, 
binary parameter measurement precisions become non-trivial 
once the Bayes factor favoring the presence of a signal ex¬ 
ceeds our threshold value of 100. Rosado et al. (2015) inves¬ 
tigated the likely properties of the first detectable continuous 
GW source in IPTA and SKA (Janssen et al. 2015) data, ob¬ 
serving that the detection probability favored massive, nearby 
binaries with orbital frequencies <10 nHz. However, the au¬ 
thors did not consider eccentricity. From our results, we see 
that a threshold detection will provide an eccentricity mea¬ 
surement precision of ^0.3 (considering the 68% credible 
region width), which may allow the eccentricity to be suffi¬ 
ciently constrained as to perform inference on plausible envi¬ 
ronmental coupling influences (such as 3-body stellar scatter¬ 
ing or circumbinary-disk interaction) which drove the binary’s 
orbital evolution. Doing so will shed light on the astrophysics 
and environment of the binary’s host galactic nucleus. 


6.3. Circular-model penalty 
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We now address the detection penalty one might incur by 
searching for an eccentric binary signal with a circular wave¬ 
form model, which was also investigated in Zhu et al. (2015) 
using frequentist methods. We firstly investigate this using 
signal injections in Type 1 data for a variety of injected binary 
frequencies and eccentricities. The matched-filtering SNR for 
a circular (monochromatic) template in data with an eccen¬ 
tric signal is computed, and compared to the optimal SNR of 
the same eccentric signal. The resulting statistic, Pdrc/Popt. 
is a measure of the effectualness of the circular template in 
representing the eccentric signal (Buonanno et al. 2009). At 
each eccentricity, the SNR is averaged over 10^ binary ori¬ 
entations and locations, and maximized over the frequency 
of a monochromatic template. The matched-filtering SNR is 
computed in three different ways: (a) as a coherent SNR for 
the entire pulsar array, maximized over the monochromatic 
template frequency; (b) as a coincident SNR, with the SNR 
in each pulsar independently maximized over the monochro¬ 
matic template frequency, and then added in quadrature to 
give the full array statistic; (c) as a coincident SNR, with the 
SNRs added in quadrature to give the full array statistic, but 
demanding a common frequency for the monochromatic tem¬ 
plate. 

Our results for Type I data are shown in Fig. 10 for or¬ 
bital frequencies beyond the region of peak PTA sensitivity 
(> 5 nHz). For cases {a) and (h), the favored monochro¬ 
matic template frequency is twice the orbital frequency until 
e ^ 0.5-0.6, and incurs an increasingly harsh SNR penalty 
as the eccentricity of the signal is increased. However, be¬ 
yond e ^ 0.5-0.6 the SNR recovers slightly, since the tem¬ 
plate frequency now favors the fundamental harmonic of the 
signal, which is lower and closer to the region of peak PTA 
sensitivity. This is seen even more clearly in case (c), where 
there is a common monochromatic template frequency across 
all pulsars when constructing the coincident SNR. The loss 
in SNR is slightly greater in {a) than in (h) and (c), since in 
the former we require signal coherence amongst all pulsars in 
the array. The behaviour found for these three cases is likely 
pessimistic, since in real matched-filtering searches the SNR 
is maximized over all template parameters rather than just the 
frequency. Lower orbital frequency SNR curves exhibit simi¬ 
larly increasing penalties as the eccentricity is raised, but the 
trends are not as smooth. 

We now investigate the circular-model penalty in terms 
of Bayesian model selection, by injecting fiducial signals of 
varying eccentricity and SNR into Type II data. The results 
are shown in Fig. 11, with the quoted Bayes factors corre¬ 
sponding to circular versus eccentric signal models. The top 
panel shows the results for the full Bayesian pipeline, whilst 
the bottom panel shows results for the eccentric statis¬ 
tic. Both techniques exhibit the same general trends: (1) at 
eccentricities ^0.1 the eccentric model receives an Occam 
penalty, resulting in the circular model being slightly favored, 
although not decisively so; (2) as the signal eccentricity in¬ 
creases so does the circular-model penalty, until the models 
are most easily discriminated at ^ ^ 0.5-0.6; (3) at higher 
eccentricities the signal model is being dominated by the fun¬ 
damental harmonic, allowing the circular model to function 
as a better approximation to the injected signal than at inter¬ 
mediate eccentricities, resulting in a reduction in the circular- 
model penalty. These trends are qualitatively similar to those 
found by Zhu et al. (2015) in Fig. (13) of their paper, however 
our eccentric search strategies exhibit superior performance at 
high eccentricity by virtue of modeling the distribution of sig- 



Figure 10. Ratio of circular-template matched-filtering SNR to optimal SNR 
for signals with various orbital frequencies and eccentricities. Case (a) shows 
results for a coherent array SNR. Case (b) shows results for a coincident ar¬ 
ray SNR with independently maximized template frequencies in each pulsar. 
Case (c) shows results for a coincident array SNR with a common monochro¬ 
matic template frequency. Further details and discussion are provided in the 
text. 



Figure 11. Bayes factors, B, for circular versus eccentric signal models when 
analyzing Type II datasets which have signals with varying injected eccen¬ 
tricity. The upper panel shows results for the full Bayesian analysis, while 
the lower panel shows results from mapping the posterior distribution of the 
intrinsic parameter space with the eccentric statistic. 

nal harmonics with in-code adaptation (see Fig. 2) rather than 
just including the lowest two harmonics. A key result of our 
analysis is that we require SNR > 8 in order for an eccentric 
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signal model to be correctly discriminated and favored when 
the true eccentricity is greater than ^0.3. 

7. CAVEATS & FUTURE DIRECTIONS 

The analysis and results presented in this paper have re¬ 
lied on several assumptions. We discuss these here, and the 
prospects for relaxing these caveats in future work. 

7.1. Prospects for including the pulsar term 

In the majority of this paper, we have ignored a full treat¬ 
ment of the pulsar term signal. Since the pulsar term is re¬ 
tarded with respect to the Earth term, it will represent the bi¬ 
nary at an earlier stage of its orbital evolution, with a larger ec¬ 
centricity and smaller orbital frequency. It is now well known 
that the pulsar term aids detection prospects for continuous 
wave sources, and is crucial in breaking degeneracies between 
the binary mass and its luminosity distance by providing ex¬ 
tra information from the binary’s evolution over the lag time 
between the Earth and pulsar term signals (Corbin & Cornish 
2010; Lee et al. 2011; Ellis 2013). 

Being able to model the orbital evolution of the binary, and 
constrain the properties of this evolution through continuous 
GW searches with PTAs, will provide a unique opportunity 
to probe the influence of other non-GW driving mechanisms. 
For example, the rate at which the binary orbital frequency, 
E, is driven by GWs, stellar scattering, and circumbinary disk 
interactions, scales as oc oc oc respectively 

(Sesana 2013). If we can include parametrized models of the 
rate of binary evolution in constructing full Earth and pul¬ 
sar term signal models in a Bayesian or frequentist search, 
then we will be able to make statements about the relative 
importance of the aforementioned mechanisms. This in it¬ 
self may provide clues as to how binaries are driven to sub¬ 
parsec orbital separations after dynamical friction in post¬ 
merger galaxies becomes inefficient, thereby adding to our 
knowledge of how the final parsec problem (Milosavljevic & 
Merritt 2003) is ameliorated. 

For now, we estimate the degree to which employing only 
the Earth term in searches is sub-optimal for detection. We 
compute the matched-filter SNR for an Earth term template 
applied to a full signal (including the pulsar term), and com¬ 
pare this to the optimal SNR for the full signal. In construct¬ 
ing the pulsar term component of the signal, we evolve the 
orbital parameters of the binary backwards in time accord¬ 
ing to Eq. (6) and the procedure outline in Sec. 4, where we 
assume all pulsars lie at a distance of 1 kpc from the Earth. 
We assume all orbital evolution is GW driven. No pericenter- 
direction evolution or orbital-plane precession is considered, 
and we do not evolve the binary during the pulsar observation 
timespan of 10 years. To ease the computational burden, we 
use a sub-array of 6 pulsars spread across the sky, averaging 
the SNR over 10^ binary locations and orientations. The ob¬ 
servational cadences and timing baselines of the sub-array are 
of Type 1 variety. 

The results are shown in Table 1 for a variety of Earth 
term orbital frequencies and eccentricities. As the orbital fre¬ 
quency and chirp mass are increased, the ratio of the Earth 
term SNR to the full SNR tends to grow with eccentricity. 
This is because higher mass, frequency, and eccentricity bi¬ 
naries are driven rapidly via GW emission, which in the most 
extreme cases leads to signals with pulsar term frequencies 
which are so far below the PTA sensitivity band that an Earth 
term template becomes an excellent approximation to the full 



Figure 12. Exclusion regions in binary eccentricity and orbital frequency as 
a function of chirp mass, corresponding to parameter combinations where 
the fundamental (dashed black lines) and second harmonic (solid black lines 
on the boundary of shaded exclusion regions) of the orbital frequency evolve 
during T = 10 years by more than the PTA frequency resolution, A/ = l/T = 
3.2 nHz, rendering the assumption of binary non-evolution invalid. 

signal. Even at fixed orbital frequency and eccentricity, the 
effect of increasing binary chirp mass is to raise the efficacy 
of an Earth term only template. However, care must be taken 
in the intermediate case, when we have moderate eccentric¬ 
ities, frequencies, and masses, which generate pulsar term 
signals that remain in the PTA band, and whose spectrum 
of GW frequencies may exceed the fundamental harmonic of 
the Earth term signal. The worst matches between signal and 
template occur for low mass, low eccentricity systems with 
orbital frequencies close to the region of peak PTA sensitvity 
(^ 5 nHz) - the combination of high array sensitivity and neg¬ 
ligible orbital evolution leads to pulsar term signals close to 
this region of peak sensitivity, and thus very poor matches 
(which are sometimes negative since we employ a coherent 
SNR). In general, the pulsar term increases the signal detec¬ 
tion prospects, but confusion may arise between different har¬ 
monics in the Earth and pulsar terms, which would harm pa¬ 
rameter estimation efforts. So long as the Earth and pulsar 
terms remain distinguishable, we will learn more about the 
system parameters from the pulsar term’s inclusion. Future 
work should study the prospects for incorporating the pulsar 
term in eccentric binary search strategies, and investigate the 
rich science that can be mined from having access to snap¬ 
shots of the binary evolution from thousands of years in its 
past. 

7.2. Binary orbital evolution during the observation 
timespan 

We now test the assumption of binary non-evolution over 
typical PTA observation timespans. For different initial Earth 
term parameter choices {fA^F^e}, we numerically evolve a 
binary forward in time by 10 years according to Eq. (6). Fig¬ 
ure 12 shows exclusion regions in parameter space where 
the fundamental and second harmonic of the orbital fre¬ 
quency evolve by more than the PTA frequency resolution, 
A/ = \/T = 3.2 nHz, which may render the approximation of 
binary non-evolution within our observing window invalid. 
The second harmonic will dominate the signal for low ec¬ 
centricities whilst the fundamental harmonic will dominate at 
higher eccentricites. 

A more rigorous way of testing this is to investigate how 
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Table 1 

The matched-filter SNR for an Earth term template is compared against the optimal full signal SNR to construct Peauh/Pfuii- At each Earth term orbital 
frequency and eccentricity, we evolve a binary backwards in time by L(1 -l-r2 • w) to construct the pulsar term waveform, where L = 1 kpc for all pulsars. 


Orbital frequency [nHz] 


e = 

0.0 



e = ( 

3.25 


Eccentricity 
e = 0.50 



e = 

0.75 



e = 

0.90 


10^ 

M [Mq] 
10* 10® 

10'® 

10^ 

M[Mq} 
10* 10® 

10'® 

10^ 

M[Mq} 
10* 10® 

10'® 

10^ 

M[Mq] 
10* 10® 

10'® 

10^ 

M [Mq] 
10* 10® 

10'® 

0.1 

0.76 

0.76 

0.76 

0.76 

0.73 

0.73 

0.73 

0.73 

0.51 

0.51 

0.51 

0.51 

0.12 

0.12 

0.12 

0.12 

0.13 

0.13 

0.13 

0.22 

0.5 

0.78 

0.78 

0.78 

0.78 

0.66 

0.66 

0.66 

0.66 

0.33 

0.33 

0.33 

0.34 

0.15 

0.15 

0.16 

0.41 

0.07 

0.07 

0.36 

0.87 

1.0 

0.63 

0.63 

0.63 

0.64 

0.47 

0.47 

0.47 

0.49 

0.25 

0.25 

0.25 

0.36 

0.10 

0.10 

0.16 

0.76 

0.05 

0.12 

0.69 

0.96 

5.0 

-0.03 

-0.03 

0.07 

0.65 

-0.02 

-0.02 

0.13 

0.72 

0.0 

0.01 

0.38 

0.87 

0.02 

0.16 

0.81 

0.99 

0.11 

0.71 

0.99 

1.0 

10.0 

-0.01 

0.02 

0.59 

0.62 

-0.03 

0.01 

0.6 

0.73 

-0.04 

0.10 

0.64 

0.93 

0.03 

0.50 

0.87 

1.0 

0.34 

0.81 

0.99 

1.0 

50.0 

0.21 

0.68 

0.63 

0.42 

0.3 

0.68 

0.56 

0.95 

0.53 

0.67 

0.81 

0.99 

0.68 

0.66 

0.99 

1.0 

0.60 

0.97 

1.0 

1.0 

100.0 

0.66 

0.65 

0.55 

0.34 

0.67 

0.64 

0.65 

0.99 

0.68 

0.52 

0.93 

1.0 

0.64 

0.81 

0.99 

1.0 

0.73 

0.99 

1.0 

1.0 


this assumption affects our ability to perform parameter esti¬ 
mation. If the non-evolution model performs well within the 
range of expected SNR, such that the systematic bias intro¬ 
duced via our assumption of binary non-evolution is smaller 
than statistical errors, then we can judge the model to be an 
excellent functioning approximation. More formally, we want 
to satisfy the indistinguishability criterion (Cutler & Vallis- 
neri 2007; Creighton & Anderson 2012): 

mt)\5s{t))<\, (40) 

where ds{t) corresponds to the difference between the ap¬ 
proximated residuals in the non-evolution model and the true 
residuals. Satisfying the inequality in Eq. (40) approximately 
corresponds to the systematic errors arising from modeling 
bias being smaller than statistical measurement errors. The 
tolerance SNR, ptoi., above which systematic errors from in¬ 
sufficient template accuracy may exceed statistical measure¬ 
ment errors, and thus become problematic, is given by 



e 


2 (S(?)|S(0) 

((5s(f)|(5s(0)’ 


(41) 


where s(0 are the true residuals (concatenated over all pul¬ 
sars) induced by a binary which may be evolving over our ob¬ 
servation timespan. To compute this, we numerically evolve 
the orbital parameters of a binary over the 10 year timing 
baseline of a Type I dataset using Eq. (6), with varying choices 
of initial orbital frequency and eccentricity. The evolved or¬ 
bit is then used to compute the pulse redshift and (via nu¬ 
merical integration) the GW-induced timing residual at each 
pulse TOA. The typical ratio of the time required to compute 
the GW signal numerically versus analytically is ^ 0(10"^), 
which is why a fully numerical approach is clearly intractable 
at present. 

The tolerance SNR is shown in Eig. 13 as a function of 
binary eccentricity, orbital frequency, and chirp mass. We 
choose a cutoff value of the tolerance SNR equal to 10 since 
this may correspond to realistic values of the SNRs of first 
PTA detections of single GW sources after ^10 years of 
IPTA and SKAl activity (Rosado et al. 2015). If our model 
can be successfully applied to real signals above this cutoff 
value, then we conclude that the treatment used in this pa¬ 
per is valid well into the era of first PTA detections. We see 
that at a binary chirp mass of IO^Mq the tolerance SNR is 
above 10 for most frequencies and eccentricities, indicating 
that the assumption of non-evolution is valid. The approxi¬ 
mation begins to break down at higher eccentricities and fre¬ 
quencies (^ 5 X 10“^ Hz) where the rate of binary evolution 
is higher. At IO^Mq our model is appropriate at all eccentric¬ 


Figure 13. The tolerance SNR, ptoi.^ for a range of binary eccentricities, or¬ 
bital frequencies, and chirp masses is shown. This indicates the SNR above 
which systematic parameter errors (which occur by keeping binary parame¬ 
ters fixed over the 10 year PTA timing baseline) may exceed statistical mea¬ 
surement errors. 


ities for frequencies lower than 10“^ Hz, however the toler¬ 
ance SNR for F = 10“^ Hz drops below cutoff at ^ ^ 0.7, and 
at higher frequencies the assumption of binary non-evolution 
is inappropriate. Einally, for the most massive binaries with 
M = IO^^Mq, the tolerance SNR remains above cutoff for or¬ 
bital frequencies lower than 5 x 10“^ Hz at all eccentricities, 
while at 5 X 10“^ Hz the tolerance SNR only drops below 10 
at ^ ^ 0.6. 

Therefore, our assumption (which has been shared by all 
other authors in this field) of binary non-evolution over typical 
PTA timing baselines is appropriate for most frequencies at 
or below the region of peak PTA sensitivity. The approxima¬ 
tion only begins to break down for the most massive systems 
above orbital frequencies of ^ 5 x 10“^ Hz and eccentrici¬ 
ties of 0.6, allowing the signal model and analysis techniques 
developed in this article to be applied to real data with robust 
outcomes. Euture studies are required to investigate faster and 
more tractable strategies for modeling the orbital evolution of 
high mass, high frequency, and high eccentricity binaries over 
PTA timing baselines. 

8. CONCLUSIONS 

PTAs are uniquely suited to explore the dynamical evolu¬ 
tion of SMBHBs before and after they decouple from their as- 
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trophysical environments to become dominated by GW emis¬ 
sion. An increasing number of studies tend to suggest that the 
mechanisms that may drive SMBHBs to small orbital separa¬ 
tions could also lead to an increase in binary eccentricity that 
will be detectable in the frequency band of PTAs. Extracting 
this information from real data will substantially increase our 
understanding of the mechanisms that lead to the formation, 
hardening and eventual coalescence of SMBHBs. In this arti¬ 
cle we have introduced several tools to address this issue. We 
have developed a robust, accurate and computationally effi¬ 
cient Bayesian pipeline to explore the feasibility of detecting 
and reconstructing the astrophysical parameters of eccentric 
SMBHBs in PTA data, and have developed for the first time 
an eccentric statistic that is, by construction, suitable to 
study systems of arbitrary eccentricity. 

We have used these tools to determine the accuracy with 
which a simulated eccentric signal could be reconstructed, 
and have conclusively shown that the recovered and injected 
parameters are completely consistent. Our prior knowledge 
of the eccentric binary parameter space will begin to be up¬ 
dated by data once the SNR of the associated binary’s GW 
signal exceeds ^ 7. We have also shown that the automated 
waveform generation algorithm, which determines the num¬ 
ber of harmonics needed to ensure that the modeled GW sig¬ 
nal reproduces the full numerical solution with an accuracy 
better than 99.999%, prevents computational inefficiencies in 
the pipeline. 

The influence of binary eccentricity on PTA single-source 
detection prospects was also considered. Assuming that the 
sensitivity peak of a PTA to continuous wave sources is lo¬ 
cated at a GW frequency /o, we have shown that eccentricity 
will enhance the detection prospects of SMBHBs with orbital 
frequencies < /q. This is because the signal spectrum of ec¬ 
centric binaries is distributed into higher harmonics of the or¬ 
bital frequency than in the case of a circular binary, leading to 
components of the signal being located in the region of maxi¬ 
mum PTA sensitivity. On the other hand, binaries with orbital 
frequencies > fo will undergo an SNR attenuation because the 
signal power is shifted to higher frequencies where the PTA 
sensitivity is poorer and dominated by TOA measurement er¬ 
rors. In summary, systems with signals which are below band 
in the circular case get pushed into band through increasing 
eccentricity, while systems that are optimally located in fre¬ 
quency for the circular case get pushed out of band by eccen¬ 
tricity. 

We found that applying a circular waveform model in the 
analysis of data with increasingly eccentric binary signals in¬ 
curs an SNR penalty which grows with eccentricity, and is 
^ 60% at worst case for coherent and coincident analyses. 
This was also investigated in a Bayesian context, where we 
found that SNRs greater than 8 are needed in order for an ec¬ 
centric signal model to be correctly discriminated and favored 
over a circular signal model when the true signal eccentricity 
is > 0.3. 

Several of the approximations used in the techniques pre¬ 
sented in this article were briefly investigated. We found that 
for very high mass, frequency and eccentricity binaries, an 
Earth term signal model performs just as well as a full sig¬ 
nal model incorporating the pulsar term, since the binary will 
have evolved so significantly that the pulsar term signal lies 
below band. Eurthermore, the possible bias from assuming 
binary non-evolution over a PTA observation time of 10 years 
was studied, and was found to be unimportant for moderately 
massive and eccentric systems in the era of first PTA detec¬ 


tions. 

There are several topics in continuous GW searches that 
should be addressed in the near future, including the need 
to assess the possible covariances involved in simultane¬ 
ous continuous GW searches and stochastic GW background 
searches. One can imagine that the reduction in low- 
frequency sensitivity associated with having fit for the pul¬ 
sar quadratic spindown parameters will be exacerbated by 
the concurrent search for a stochastic GW background signal 
dominated by low-frequency power. Also, the possibility of 
having multiple resolvable continuous-wave sources may lead 
to difficulties in isolating each source during the Bayesian 
searches, resulting in interesting covariances. However there 
are ongoing efforts to resolve this issue (Ellis 2015). Eurther¬ 
more, as new data is added to each PTA and combined to form 
IPTA datasets, the prospects for continuous GW source detec¬ 
tion grow stronger. Current pipelines should be tested against 
signals injected into near-future type datasets as a means to 
inform new advances in analysis procedures. Einally, the effi¬ 
cacy of performing continuous GW searches on datasets hav¬ 
ing signals composed of realistic GW source populations must 
be addressed in the near future. These areas of future study 
have not yet been investigated with circular-binary GW signal 
models, however our development in this paper of more com¬ 
plete signals models which include eccentricity will endow 
these studies with greater verisimilitude. 

Huerta et al. (2015) and this article have provided a solid 
foundation to explore in a consistent way the influence of ec¬ 
centricity on the detection and parameter estimation of SMB¬ 
HBs with PTAs. The tools presented in this article can be 
readily incorporated into all present and planned analysis 
pipelines. The toolkit introduced in these articles could be ex¬ 
tended to explore in detail what constraints may be placed on 
the various astrophysical mechanisms that can drive the dy¬ 
namical evolution of SMBHBs prior to becoming dominated 
by GW emission. 
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